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Abstract. We study atomic Josephson junctions (AJJs) with one and two bosonic 
species confined by a double-well potential. Proceeding from the second quantized 
Hamiltonian, we show that it is possible to describe the zero-temperature AJJs 
microscopic dynamics by means of extended Bose-Hubbard (EBH) models, which 
include usually-neglected nonlinear terms. Within the mean-field approximation, the 
Heisenbcrg equations derived from such two-mode models provide a description of 
AJJs macroscopic dynamics in terms of ordinary differential equations (ODEs). We 
discuss the possibility to distinguish the Rabi, Josephson, and Fock regimes, in terms of 
the macroscopic parameters which appear in the EBH Hamiltonians and, then, in the 
ODEs. We compare the predictions for the relative populations of the Bose gases atoms 
in the two wells obtained from the numerical solutions of the two-mode ODEs, with 
those deriving from the direct numerical integration of the Gross-Pitaevskii equations 
(GPEs). Our investigations shows that the nonlinear terms of the ODEs are crucial 
to achieve a good agreement between ODEs and GPEs approaches, and in particular 
to give quantitative predictions of the self-trapping regime. 
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1. Introduction 

The prediction [1] of Bose-Einstein condensation (BEC) and the experimental 
achievement of BEC [2] has played a crucial role for theoretical and experimental 
developments in the physics of ultracold atoms. The study of the atomic counterpart [3, 
4, 5, 6, 7] of the Josephson effect which occurs in superconductor-oxide-superconductor 
junctions [8] - which is an example of macroscopic quantum coherence - represents one 
of these developments. Albiez et al. [9] have provided the first experimental realization 
of the atomic Josephson junction (AJJ) previously analyzed theoretically in a certain 
number of papers [3, 4, 5, 6, 7]. In 2007 Gati et al. [10] reviewed the experiment by 
Albiez et al. [9] and compared the experimental data with the predictions of a many- 
body two-mode model [11] and a mean-field description. In the above references the 
analysis of AJJs physics is carried out in the presence of a single bosonic component. The 
possibility to tune intra- and inter-species interactions [12, 13] by means of the Feshbach 
resonance technique makes possible to study of AJJs with two bosonic species trapped 
together by double-well potentials and to use BECs mixtures as powerful instruments 
to investigate quantum coherence and nonlinear phenomena, with particular attention 
to the existence of self-trapped modes and intrinsically localized states. 

In the superfluid regime the dynamics of the relative populations and relative phases 
of the Bose condensed atoms can be described by Josepshon's two-mode equations, which 
are ordinary differential equations (ODEs), see for example Refs. [5, 14, 15, 16, 17]. 
This description is achieved in the presence of a confining double-well potential, with 
a single bosonic component [5] and also with bosonic mixtures [14, 15, 16, 17, 18]. 
One of the most interesting aspects of AJJs analysis is to compare the predictions 
deriving from the ODEs with the ones obtained from the Gross-Pitaevskii equations 
(GPEs). For single component condensates, Salasnich et al. [6] have shown that a 
good agreement exists between the results obtained from the GPE and those of the 
ODEs. Similar agreement was obtained in [7] for AJJs realized with weakly interacting 
solitons localized in two adjacent wells of an optical lattice. However, the situation may 
be quite different for multicomponent condensates, due to the interplay of intra- and 
inter-species interactions which enlarges the number of achievable states (for instance, 
mixed symmetry states can exist only in the presence of the inter-species interaction) as 
well as their stability, giving to the system many more dynamical possibilities. Recently 
it was shown that for the two components case the integration of the ODEs allows to 
predict the analogous of the macroscopic quantum self-trapping phenomenon observed 
in AJJs with one bosonic component [15, 17]. This phenomenon has been discussed for 
a two-components nonlinear Schrodinger model with a double-well potential by Wang 
and co-workers [19]. More recently, a comparison between the reduced ODEs system 
and the full GPE dynamics was performed, showing that, for various conditions, a good 
agreement exists between the two kinds of predictions [17]. 

The aim of the present work is to analyze how the accuracy of the two-mode 
approximation can be improved by taking into account the usually-neglected nonlinear 
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terms. These terms derive from the overlaps between wave functions localized in 
different wells. Both for single component and for two components AJJs - introduced 
in the second section - we proceed from a full second quantized description of the 
system. In Sec. Ill we describe the system by the extended Bose-Hubbard (EBH) 
Hamiltonian. In the single component case, the EBH Hamiltonian is the two-sites 
restriction of the Hamiltonian considered in Refs. [20, 21] to analyze bosons loaded in 
one dimensional optical lattices. In the two species case, the EBH Hamiltonian is the 
extended version of the one considered by Kuklov and Svistunov in Ref. [22] to study 
the counterfiow superfluidity of two-species ultracold atoms. We note that the study 
of the two components bosonic system proceeding from a pure quantum approach is a 
subject of wide interest. In fact, this topic is dealt with in certain regions of the phase 
space in Ref. [23] and in the case of hardcore bosons as discussed in Ref. [24] . 

The EBH Hamiltonian sustains the dynamics of the single-particle operators via the 
Heisenberg equations of motion [25, 26]. By performing the mean-field approximation 
on the single-particle operators of each component, the improved ODEs are achieved. 
In the third section we also discuss how it is possible to distinguish the Rabi regime, 
the Josephson regime, and the Fock regime. This analysis is carried out in terms of the 
macroscopic parameters involved in the EBH Hamiltonians and, then, at the right hand 
sides of the improved ODEs as discussed for single AJJs in Ref. [27]. In Sec. IV we 
write down the GPEs for the one and the two components AJJs. Here we compare the 
results obtained by numerically integrating the GPEs with the predictions obtained by 
numerically solving the improved ODEs. Moreover, in the fourth section we plot the 
phase-plane portraits of the dynamical variables fractional imbalance-relative phase. 
Finally, in Sec. V we draw our conclusions. 

2. The system 

We consider two interacting dilute and ultracold Bose gases denoted below by 1 and 2. 
We suppose that the two gases are confined in a double-well trap produced, for example, 
by a far off-resonance laser barrier that separates each trapped condensate in two parts, 
L (left) and R (right). We assume, moreover, that the two condensates interact with 
each other and that the trapping potential V trap (r) for both components is taken to be 
the superposition of a strong harmonic confinement in the radial (x-y) plane and of a 
double-well (DW) potential in the axial (z) direction. We model the trapping potential 
as: 

V trap (r) = ^(x 2 + y 2 ) + V DW (z) , (1) 

where is the mass of the ith component. For simplicity we take uj\ = oj<i = oj. For 
symmetric configurations in the z direction, we take - for the ith species - the double- well 
in Eq. (1) as 

V D w(z) = V L (z) + V R (z) 
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that is the combination of two Poschl- Teller (PT) potentials, Vl(z) and Vr(z), centered 
at the points — zq and Zo, and separated by a potential barrier which may be changed 
by varying b (see Fig. 1). We use PT potentials only for the benefit of improving 
accuracy in our numerical GPEs calculations (see the fourth section), taking advantage 
of the integrability of the underlying linear system. We remark, however, that our 
results apply to a generic double-well potential. Eigenvalues and eigenfuctions of the 
Poschl- Teller potential for a single well are known analytically. The wave functions of 
the ground state of V a (z) (a = L,R), centered around — z (+Zq) are [28] : 
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The constant A, in Eq. (3), ensures the normalization of the wave function in each well. 
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Figure 1. The double- well potential (2) as a function of z for zq — 3 and different 
values of b. The dot-dashed line corresponds to b = 0.7, the continuous line corresponds 
to b = 1, and the dashed line corresponds to b = 1.3. Lengths are measured in units 

of a±_ i = \ and energies in units of hw . 
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3. The second quantization Hamiltonian 

To describe our system at zero-temperature, we proceed from the second quantized 
Hamiltonian, which reads 

,2 



H=J2 /rf 3 r^(r)(-^-V 2 + Wr))^W 

1=1,2"^ * 

+ E | / dh^rMv)h(r)Mr) 

i=l,2 J 

+ g 12 J dW^r^r^r)^!-) , (4) 

where V trap (r) is the potential (1). The coupling constants </j and g\ 2 are the intra- and 
inter-species atom-atom interaction strengths, respectively. These constants are given 
by 

4:7rh 2 ai 

9i = , ( 5 ) 

rrii 

9i2 = , 6 

m r 

where the reduced mass m r is equal to m 1 m 2 /(m 1 + m 2 ). Eqs. (5) and (6) relate 
the two coupling constants to the respective s-wave scattering lengths, a« and a\ 2 . 
In the following, we shall consider both and g\ 2 as free parameters, due to the 
possibility of changing the s-wave scattering lengths Qj and a 12 by the technique of 
Feshbach resonances. In the following, we will neglect the mass difference between the 
two bosonic components of the mixture, as for example in Ref. [13], and assume that 
mi = m 2 = m. In Eq. (4), the field ^(r) (\&|(r)) destroys (creates) a boson of the iih 
species at the point r, and obeys the usual bosonic commutation relations. We expand 
the field operator ^(r) in terms of operators a aji (a' ai ) - destroying (creating) a boson 
of the zth species in the well a = L,R - according to: 

*<(r) = E ®*A r Xi , ( 7 ) 

a=L,R 

where a's and a^'s satisfy the usual boson commutation relations and the functions $ Qj j 
form an ortho normal set. Due to the form (1) of the trapping potential, 3> aj j(r) can be 
decomposed as 

$ a ,i(r) = Wi(x)wi(y)(j) at i(z) , (8) 

where Wi{x) and Wi(y) are the ground state wave functions of the harmonic oscillator 
potentials miU 2 x 2 /2 and m,iU! 2 y 2 /2, respectively. The functions 4>L,i(z) and <f>R t i(z) at 
right hand side of Eq. (8) are two functions well localized in the left and right well, 
respectively. These functions are real and orthonormal. The functions 4>L,i{z) and 
<l>R,i{ z ) can be determined following the same perturbative approach as in Ref. [17]. 



Nonlinear quantum model for atomic Josephson junctions with one and two bosonic speciesG 



Under the same conditions, these functions may be written in terms of the 4>(l,i,pt){ z ) 
and 4>{R,i,PT)( z ) of Eq. (3) as: 




where s = f_™ dz (j>(L,i,PT)(z)<j>(R,i,pr)(z)- 
3.1. AJJs with a single bosonic species 

Let us start our analysis by considering the presence of a single bosonic component. 
In this case, the inter-species coupling constant (6) is equal to zero. We use the field 
operator expansion (8) in the second quantized Hamiltonian (4). The AJJs microscopic 
dynamics is controlled by the EBH Hamiltonian [20, 25, 26]. The EBH model, by 
omitting the species index i, is described by the Hamiltonian 



Here the number of particles in the ath well. E a are the energies of the 

two wells, U a > are the boson-boson repulsive interaction amplitudes, and K is the 
tunnel matrix element, which is the Rabi oscillation energy in the case of a model with 
U a equal to zero. The parameter K c is the induced collisionally hopping amplitude, V 
is the density- density bosonic interaction amplitude, and K p describes the pair bosonic 
hopping [20]. By using the decomposition (8) and the explicit form of w(x) and w(y), the 
macroscopic parameters (10) may be shown to be related to the intra-species coupling 
constant (5) and to the other microscopic parameters (the mass and the frequency of 
the harmonic trap) by the formulas 
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K v = \, (11) 

where g = — ^-y. We observe that the first two lines of Hamiltonian (10) involve only 

the overlaps between a 's localized in the same well, see [5]. The third and fourth lines of 
the Hamiltonian (10) include also the overlaps between a 's localized in different wells, 
see [27]. Proceeding from the Hamiltonian (10), we write down criteria to individuate 
different oscillations regimes sustained by the AJJs dynamics. To this end, as discussed 
in Refs. [4, 27], we express the Hamiltonian (10) in terms of the following operators: 
1 

Jx = 2^L^L ~ a R a R) 

J* = \{a ] L a R + a R a L ) , (12) 

and the SU(2) algebra invariant J 2 = (N/2)(N /2 + 1), with N being equal to hi + tir 
[29]. 

We assume that the two potential wells are symmetric, E\ = E R = E and 
Ul = Ur = U . Neglecting constant terms, and using the fact that N ^> 1, we get 

H = (U - V) J 2 - 2 (K - K c N) J z + VJl . (13) 

Here, if the condition (U — V) ^> V is verified, we can consider only the terms in 
and J z . Then by defining the parameter R as 

R= \ U - V) \ (14) 
(K — K C N) 1 ; 

we are able (see Ref. [4]) to distinguish the three following regimes 

• Rabi: R < 1; 

• Josephson: 1 < R < N 2 ; 

• Fock: N 2 < R. 

In the Rabi regime the bosons are in a coherent state and oscillate with a frequency 
given simply by the energy difference between the ground state and the first excited 
state associated to the double-well potential. In the Josephson regime the bosons are 
in a coherent state and oscillate with a frequency which depends on the parameters U, 
K c , and V. Moreover, if the interaction strength is sufficiently large the self-trapping 
takes place. In the Fock regime the bosons are in a Fock state characterized by the 
suppression of number fluctuations. Now, we observe that the Hamiltonian (10) can 
be viewed as the two sites restriction of the Hamiltonian considered in Refs. [20, 21] 
within the study of bosons loaded in one dimensional optical lattices. In particular, in 
Ref. [20] it is shown that within the Fock regime two regions open up. To this end we 
denote by A the energetic gap between the Fock state with iV bosons per well and the 
Fock state with (N + 1) per well [20] 

A = 2(E + U N ) . (15) 
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Then, when 

|4A-(2JV + 1)V| > V (16) 

we have a pure Mott insulating phase (PMI), driven by the density-density on-site 
interaction. When 

|4A- (2iVo + 1) V\ < V (17) 

we have a Density-Wave Mott insulating (DWMI) regime, driven by the nearest- 
neighbors interaction [20]. Note that the DWMI phase is characterized by number 
fluctuations suppression as well. 

At this point, we remark that we are interested in determining the fully coherent 
dynamical oscillations of population of the Bose condensed atoms between the left and 
right wells. Then, we proceed from the Heisenberg equations of motion for the model 
Hamiltonian (10). These equations of motion control the temporal evolution of a a . We 
observe that in the superfluid regime the system is in a coherent state and the following 
mean-field approximation [25] 

(a a ) = \/N~e^g{i9 a ) 

«) = N a (18) 

can be performed. The averages involved in Eq. (18) are evaluated with respect to the 
coherent state. Under the assumption of symmetric wells and by inserting the mean- 
field approximation (18) into the aforementioned Heisenberg equations of motion, we 
get 

2(K-K C N) , 

z(t) = — ^ c 1 y/l-z 2 (t) Bin 0(t) 

+ ^(l-z 2 (t))sin20(t) 

• = U-V 2(K-K C N) McceO® 
h h y/l-z 2 (t) 

VN 

-^-z(t)cos26(t) , (19) 

where N = Nl+Nr is the total number of bosons, and z = (Nl—Nr)/N and 9 = 9r — 9l 
are, respectively, the fractional imbalance and the relative phase. 

3.2. AJJs with two bosonic species 

In this subsection we shall consider AJJs in the presence of two interacting bosonic 
components. In this case both the coupling constants (5) and (6) are finite, and the two 
mode EBH model is described by the Hamiltonian 

H = ^2 H(EBH,i) + H12 . (20) 
i=l,2 

The Hamiltonian H(ebh,i) is the single component Hamiltonian (10) written in terms 
of the operators a a ,i and aJ ai . The parameters E®, U a , K, K c , V, K p , and the function 
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(p a will read E® { , U a> i, K iy K Cji , Vj, K p>i , and 4>a,ii respectively. The microscopic 
quantities referred to a single bosonic component will be modified according to the 
same prescription. Under the hypothesis of symmetric wells, the coupling Hamiltonian 
H12 reads: 

H\2 = Ui2(d^ L1 a) L2 aL,idL t 2 + ^r^r^R^R?) 

+ Vui^L^R^L^aR^ + a ^L,2®R,l a L,2(iR,l) 

+ K P 

,l2\ a L i a L 2 a ^,2 a -R,l + a R i a R 2 CL L,20'L,1 
+ ^L,l^R,2^L,2dR^ + Ojr^l^R^L,!) 

+ K Ci i 2 (a/ L1 hL,2aR,i + d Rjl h Li2 d Lil 

+ a , L 2^ l L ) \0'R,2 + ^R^L,ldL, 2 

+ a) L2 h R ^aR2 + d R2 fiRidL t 2 

+ 0^71^20^,1 + a R1 nR !2 a Ltl ) . (21) 

In Eq. (21), U12 is the inter-species interaction amplitude between bosons localized 
in the same well, and V12 is the inter-species interaction amplitude between bosons 
localized in different wells. The quantity i^ p ,i2 is the inter-species pair hopping (hopping 
of particle-particle or hole-hole pair made up of bosons of different species); i^ Cj i2 is the 
amplitude of the inter-species collisionally induced hopping. By using the decomposition 
(8) and the explicit form of Wi(x) and Wi(y), the aforementioned parameters are shown 
to be related to the inter-species coupling constant (6) by: 



Ui2 = gi2 dz (<f) a j(z)y 
J —00 
r+00 

V12 = 9l2 / dz (0 Qji (2)) 2 (0^j(2)) 2 



00 

' 00 



K c ,i2 = gu / dz {(j) a ^{z)f((t)pj{z)) 
J —00 

= V12 , (22) 

9l2 

where g i2 = — r-s 5 — -. Note that we are considering both the overlaps between 

+ «i, 2 ) 

4>aS localized in the same well (JJ a ,i and U12) - that are the only terms taken into 
account in Ref. [17] - and the overlaps between a 's localized in different wells (V12, 
Kp,i2, Kc,n)- We observe that, in general, due to the presence of the parameters (22) 
the identification of different oscillation regimes proceeding from the Hamiltonian (20) 
is not immediate as for single component AJJs. Nevertheless, under certain conditions 
we are able to write down criteria to select the different regimes sustained by the two 
components AJJs dynamics. First, let us focus on the case in which only the overlaps 
between 0q,'s localized in the same well are considered. If certain relations exist between 
the intra- and the inter-species interactions amplitudes, we can recognize the two-species 
corresponding of the Rabi, Josephson and Fock regimes discussed in the case of single 
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component AJJs. For each component i, we define the quantity 73 as 

UiNi 

H = —ft- ■ (23) 
We recognize the following "weak-coupled" Rabi, Josephson, and Fock regimes 

• Rabi: 7, < 1, \U 12 \ ~ Uf, 

• Josephson: f <C 7^ <C iV 2 , |[/" 12 | < t^; 

• Fock: JV? < 74. 

In the Josephson regime, even if the intra-species interaction is not strong enough to 
ensure self-trapping by itself, self-trapping occurs when the inter-species interaction 
strength exceeds a crossover value. In the Fock regime the net number of atoms in the 
transport is suppressed. However, with repulsive inter-species interaction the so-called 
counterflow survives [22]. This means that the currents of the two species are equal 
in absolute values and are in opposite directions. This conductive regime is named 
super (counter)fluid phase (SCF). As discussed in Ref. [22], the system supports the 
SCF phase of the two components when 

U x + U 2 - 2 U 12 > 1 . (24) 

When the condition 

U x + U 2 = 2 U 12 (25) 

is met, a phase separation (PS) is observed in the system and the system can be viewed 
as composed by two totally independent Bose gases confined in the double-well potential. 
On a physical level, this phase separation means that one bosonic component will occupy 
the left well and the other the right well. If the inter-species interaction is attractive 
and the hypothesis Ni — N 2 = N is verified, then, when 

E/i + Z7 2 -2|E/i2|»l (26) 

a superfluid phase, in which the superfluid consists of pairs of bosons, is supported by 
the system. This phase is named superfluid paired phase [31]. 

So far we have neglected the role played by the terms deriving from the overlaps 
between Q 's localized in different wells. The presence of these terms makes the scenario 
more complicated. However, also in this situation, under certain conditions, it is possible 
to achieve a classification of the oscillations regimes. To this end, as discussed for the 
single component case, we express the Hamiltonian (20) in terms of the operators J x ^ 
J y ,i, Jz,i defined in Eq. (12) and the SU(2) algebra invariant J? = (A r i /2)(A r i /2 + 1), 
with Ni being equal to fiL,i + n Rt i. Since we are assuming symmetric potential wells, 
we can write that E Q Li = E° Ri = Ul^ = Ur^ = U,i. Neglecting constant terms, and 
using the fact that Ni 3> 1, we get 

H=(Ui- Vi)Jli — 2{K — K C)l N, - K c , 12 N 3 )J Z , 
+ Vj 2 ^ + 4(([/ 12 - V 12 )J Xjl J x ,2 + V 12 J Z:1 J Z , 2 ) 
+ U 12 (n LA n Rt2 + n Rt in Li2 ) + V\ 2 {h L ±h Ly2 + n RA h Rt2 ) . 

(27) 
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Again, if (C/j — Vj) 3> Vj, Vi 2 , and (Ui 2 — V\ 2 ) 3> V$, Vi 2 , we can consider only the terms 
in J zi , and J Xi iJ x>2 - We will assume also that Nx — N 2 = N, U\ — U 2 = t/, 
Vi = V2 = V, and iT C) i = i^ C 2 = iC C j and that the initial conditions are the same for 
both the components. In analogy to the case of a single component AJJ, we define the 
parameter R as 

_ {(U-V)+4(U vl -V 12 ))N 

(K-(K c + K cA2 )N) ■ [m) 

Again, we are able to distinguish the three regimes: 

• Rabi: R < 1; 

• Josephson: 1 < R < iV 2 ; 

• Fock: iV 2 < 

At this point, we remark that we are interested in determining the fully coherent 
dynamical oscillations of population of the two bosonic components between the left 
and right wells. Then, we proceed from the Heinseberg equations of motion for the 
model Hamiltonian (20). These equations of motion control the temporal evolution of 
a a> i. Again, by inserting the mean-field approximation valid in the superfluid regime - 
(o>a,i) = \/Naj exp(i6 l Qj j), (n at i) = N a ^ - into the aforementioned Heisenberg equations 
of motion, one gets the coupled differential equations for the fractional imbalance 
Zi = (Ni,i — Nn t i)/Ni and relative phase 6>j = 9r^ — 9l,z of the two species: 

2{K t - K c4 Ni] 



m = - h v 1 - ™ sin ^ (t) 

+ ^(l-^ 2 (t))sin20 i (t) 



2 



--^(Vnyjl-zftt) cob 6 jit) 



+ K^ 2 )N^l-zf{t) sin^(t) 

- M Zl (t) cos 2^) + Ul2 ~ Vl2 N jZj {t) 
2 



'-(V 12 ^l-z](t) cob 9 j(t) 
+ K c , 12 )Nj 



Zj(t) cos 9j(t) 



(29) 



4. Gross-Pitaevskii equations predictions: comparison with ordinary 
differential equations results 

So far we have discussed how AJJs dynamics can be described by means of the ODEs, 
i.e. Eqs. (19) and (29). We know that AJJs dynamics can be analyzed, in the mean-field 
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approximation, in terms of partial differential equations, i.e. the GPEs. This description 
can be achieved proceeding from the Heisenberg motion equations for the field operators 
\&i(r, t), (i = 1,2), associated to the Hamiltonian (4), that is 

mifi = [%,H]. (30) 

The average - denoted by (...) - of both sides of Eq. (30) evaluated with respect to the 
coherent state, provides the two coupled GPEs 
/Mr. h 2 

ifr-gf = V 2 % + \V iTap (r)+g i \%\ 2 + g ij \^j\ 2 ]% . (31) 

The macroscopic wave functions \I/i(r, t) = (^(r, £)) of interacting BECs in the trapping 
potential Vj rap (r) at zero-temperature satisfy Eq. (31). The wave function \l/,(r,t) is 
subject to the normalization condition 

y rf 3 r |*,(r,t)| 2 = ^ . (32) 

We are interested to study the dynamical oscillations of the populations of each 
condensate between the left and right wells when the the barrier is large enough so 
that the link is weak. To exploit the strong harmonic confinement in the (x-y) plane 
and get the effective one-dimensional (ID) equations describing the dynamics in the z 
directions, we write the Lagrangian associated to the GPE equations in (31) 

,2 



/ rf 3 r ([ E *, ( 4 + A^)*, 

J t=l,2 



9i 



9ij\Vi\ 2 \V 



3 I 



(33) 



where denotes the complex conjugate of and % ^ j; then, by following the 
decomposition (8) and the Gaussian approximation for the radial part of wave function, 
we adopt the ansatz 

Vi(x,y,z,t) = —= exp - X V fi{z,t) , (34) 



2a 



where the field fi(z, t) obey to f_™ dz\fi(z) | 2 = Ni, so that the normalization condition 
given by Eq. (32) is satisfied. Note that the Gaussian ansatz with the transverse width 
simply given by a±^ is reliable under very strong transverse confinements, namely when 
9i\fi\ 2 ^ 2 f\Wi [32]. By inserting the ansatz (34) in Eq. (33) and performing the 
integration in the radial plane, we obtain the effective ID Lagrangian for the field 
fi(z,t). Such an effective ID Lagrangian reads 

h 2 d 2 

~)Ji 



J i=l,2 



2irii dz 2 ' 

w*))i/<r-fi/<i 4 ] - hifiWji 2 ) 



(35) 
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Figure 2. Fractional imbalance z{t) vs. time for single component atomic Josephson 
junctions. The parameters of the double-well potential (2) are chosen to be b = 1 
and zo = 3. The dashed line represents data from the integration of GPE (37), the 
continuous line represents data from the integration of ODEs (19), and the dot-dashed 
one represents data from the integration ODEs (19) with K c = V = 0. We have set 
N = 200 and K = 4.955 x 10~ 3 . We have used the initial conditions z(0) = 0.6 and 
0(0) = 0. In the top panels (from left to right): U = 0.05 K, K c = -1.842 x 10~ 6 , 
V = 2.268 x 10~ 7 ; U = 0.1 K, K c = -3.684 x 10~ 6 , V = 4.535 x 10~ 7 . In the bottom 
panels (from left to right): U = 0.2 K, K c = -7.368 x 10~ 6 , V = 9.070 x 10~ 7 ; 
U = 0.5 K, K c = -1.842 x 10~ 5 , V = 2.268 x 10~ 6 . Time is measured in units of a;" 1 
and energies are measured in units of fvjj. 

h 2 miufa 2 ±i 

where q is given by q = 5 1 -. By varying L with respect to fa, we obtain 

2m i a z L i 2 

the ID GPE for the field ft 

8f- ft 2 8 2 f 

= + ^ + Vbwr(z) + ~ 9i \fi\ 2 + hlfiWi ■ (36) 
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Figure 3. Phase diagrams of the fractional imbalance z(t) vs. macroscopic phase 6{t) 
for single component atomic Josephson junctions. The parameters of the double-well 
potential (2) are the same as in Fig. 2. In both the panels we have set N = 200 and 
K = 4.955 x 10~ 3 . 

Left panel: the dashed line represents data from the ODE (19) with U = 0.05 if and 
K c = V = 0; the continuous line represents data from the ODE (19) with U = 0.05 K, 
K c = -1.842 x 10~ 6 and V = 2.268 x 10" 7 . 

The right panel shows the phase diagram for the self-trapping. In this panel the 
dashed line represents data from the ODE (19) with U = 0.2 K and K c = V = 0; the 
continuous line represents data from the ODE (19) with U = 0.2 K, K c = — 7.368 x 10 -6 
and V = 9.070 x 10~ 7 . Initial conditions are the same as in Fig. 2. Time is measured 
in units of and energies are measured in units of hu>. 



In the presence of a single bosonic component, g\ 2 = 0; then, the two coupled ID GPEs 
Eq. (36), omitting the species index i, reduce to 

tfi l=-^£ + i £+,/ -w + * |/|2 i / - (37 > 

Now, we observe that it is possible to write the fields fa, (i = 1, 2), by using the two-mode 
approximation as done, for example, in Ref. [17] 

fi{z,t) = $Lt{p)4>L,i{z) +i>R,i{t)<l>R > i{ Z ) 

(*) = exp (i0a,i (t)), (38) 

with 4> a ,i{.z) constructed as discussed in Sec. Ill, see Eq. (9). Then, one takes into 
account the overlaps both between </> a 's localized in the same well and between c/> Q 's 
localized in different wells. By following the same path as in Ref. [17], when the inter- 
species coupling constant g\ 2 is finite, it is possible to recover Eqs. (29) for binary AJJs, 
while for g 12 equal to zero one gets back the Eqs. (19) for single component AJJs . 

At this point - both for single component and for two components AJJs - we may 
compare the predictions of the ODEs, Eqs. (19) and (29), and those of the GPEs, 
Eqs. (37) and (36). The results of this analysis are reported in Fig. 2 for the single 
component case, and in Figs. 4, 6 for the two components case. In obtaining Fig. 2 
we have fixed the parameters b and z of the double-well potential (2). Then, by using 
the functions (9) into the third of Eqs. (11), we have obtained the tunneling amplitude 
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Figure 4. Fractional imbalance Zi(t) of the two bosonic species vs. time. The 
parameters of the double-well potential (2) are chosen to be b = 1 and zo = 3. Here, 
the dashed line represents data from the integration of GPEs (36), the continuous line 
represents data from the integration of ODEs (29), and the dot-dashed line represents 
data from the integration ODEs (29) with K Cti = Vi = K c .i2 = V\2 = (i.e. the ODEs 
analyzed in Ref. [17]). We have fixed N t = 200 and N 2 = 100. Moreover, A'j = 
K 2 = K = 4.955 x 10~ 3 , U x = U 2 = U = 0.1 K, K c>1 = K c . 2 = K c = -3.684 x lO" 6 , 
Vi = V 2 = V = 2.268 x 10" 7 . We used the initial conditions zi(0) = 0.5 = -z 2 (0) 
and 6>i(0) = 6» 2 (0) = 0. In the top panels we set U 12 = -U/20, K c<12 = -K c /20, 
V12 = -V/40, in the middle panels U 12 = -U/2, K c . 12 = -K c /2, V 12 = -V/A, and in 
the bottom panels U\ 2 = —U, K c ,\ 2 = — K c , V\ 2 = —V/2. Time is measured in units 
of (wi) _1 = (u^) -1 = uj^ 1 and energies are measured in units of hui. 
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Figure 5. Phase diagrams of the fractional imbalance Zi(t) vs. macroscopic phase 
9i(t) of the two bosonic species. The parameters of the double- well potential (2) 
are the same as in Fig. 4. In both the panels we have set Ni = 200, N 2 = 100, 
K x = K 2 = K = 4.955 x 10~ 3 , Ux = U 2 = U = 0.1K. In both the panels, the 
dashed line represents data from the ODEs (29) with Ui 2 = —U/2 and K c ,i = Vi = 
K c> i2 = V12 = 0, the continuous line represents data from the ODEs with U12 = —U /2, 
Kc,i = K c , 2 = K c = -3.684 x 10~ 6 , V x = V 2 = V = 2.268 x 10" 7 , K e ,u = -K c /2, 
V12 = —V/4. Initial conditions are the same as in Fig. 4. Time is measured in units 
of = = and energies are measured in units of fko. 
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Figure 6. Fractional imbalance Zi(t) of the two bosonic species vs. time. In the 
double-well potential (2) we set 6=1 and z = 3. In this figure the dashed line 
represents data from the integration of GPEs (36), the continuous line represents data 
from the integration of ODEs (29), and the dot-dashed line represents data from the 
integration ODEs (29) with K Ci i = Vi = K c ,i2 = V12 = (i.e. the ODEs analyzed 
in Rcf. [17]). We have fixed N x = 200 and N 2 = 100. Moreover, K 1 = K 2 = 
K = 4.955 x 10~ 3 , Ux = U 2 = U = 0.1K, K cA = K c , 2 = K c = -3.684 x lO" 6 , 
Vi = V 2 = V = 2.268 x 10~ 7 , Ui 2 = -2 U, K cA2 = -2 K c V 12 = -V. Initial conditions 
are the same as in Fig. 4. Time is measured in units of (wi) _1 = (t^) -1 = uj^ 1 and 
energies are measured in units of hxi. 
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Figure 7. Phase diagrams of the fractional imbalance Zi(t) vs. macroscopic phase 
9i(t) of the two bosonic species for the self-trapping. The parameters of the double- well 
potential (2) are the same as in Fig. 4. In both the panels we have set N\ = 200, 
N 2 = 100, K x = K 2 = K = 4.955 x 10~ 3 , U x = U 2 = U = O.IK. In both the 
panels, the dashed line represents data from the ODEs (29) with Ui 2 = —2 U and 
K C: i = Vi = K c .i2 = V12 = 0, the continuous line represents data from the ODEs 
with U12 = -2 U, K cA = K Ct2 = K c = -3.684 x 10" 6 , V 1 = V 2 = V = 2.268 x 10~ 7 , 
A' c i2 = —2K C , V12 = —V. Initial conditions are the same as in Fig. 4. Time is 
measured in units of = (^2) _1 = w^ 1 and energies are measured in units of fouj. 

K. We have keeped fixed K and we have plotted the predictions of the ODEs (19) for 
z(t) in correspondence to different intra-species interactions both when K c and V are 
zero - dot-dashed lines - and in the presence of K c and V - continuous lines; the dashed 
lines represent z(t) obtained by numerically integrating the GPE (37). In Figs. 4, 6 
we have fixed the tunneling amplitude Ki - as done previously in the single component 
case - and the intra-species interaction Ui, and we have plotted the predictions of the 
ODEs (29) for Zi(t) in correspondence to different inter-species interactions both when 
K c>i , Vi, K Cjl2 , V12 are all equal to zero - dot-dashed lines - and in the presence of 
K Ct i, Vi, K C) i2, V\2 - continuous lines; again, the dashed lines represent Zi(t) obtained 
by numerically integrating the GPEs (36). In the two top panels of Fig. 2 and in all 
the panels of Fig. 4, we have plotted the temporal evolution of the bosonic fractional 
imbalances z when they oscillate around a zero time-averaged value, i.e. (z(t)) = 0. We 
see that the usually-neglected nonlinear terms play a crucial role in order to improve 
the agreement between the GPEs the ODEs predictions. In fact, neglecting these terms, 
the solutions of ODEs and GPEs diverge rather rapidly, as shown by dot-dashed lines in 
Fig. 2 - single component AJJs - and by dot-dashed lines in Fig. 4, for two components 
AJJs. The two bottom panels of Fig. 2 show the results of our analysis when the intra- 
species interaction amplitude U is sufficiently large to induce oscillations of z(t) around 
a non zero time-averaged value, that is the self-trapping. We see that the inclusion 
within the description of the system of the usually-neglected nonlinear terms produces 
an improvement in the agreement between the ODEs and GPE predictions. In the 
two components case, from Fig. 4 we can see that the nonlinearity associated to the 
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intra-species interaction is not strong enough to induce oscillations of Zi around a non 
zero time-averaged value. Nevertheless, if the inter-species interaction is sufficiently 
large, oscillations of Zi around (Zi(t)) 7^ are observed. We have reported this kind of 
behavior for both the components in Fig. 6. From this figure we can see that, especially 
in the case of large inter-species interaction, the role played by the parameters describing 
the overlaps between Q 's localized in different wells becomes essential to improve the 
agreement between the ODEs and the GPEs predictions. Moreover, in Fig. 3 - single 
component AJJs - and in Fig. 5 and Fig. 7 - two components AJJs - we show the 
phase-plain portraits of the dynamical variables Zi and 0, for different values of the 
macroscopic parameters (11) and (22) (see Figs. 3, 5, 7 for the details). These figures 
show the comparison between the trajectories in the phases space obtained by integrating 
the ODEs in the absence of the usually-neglected nonlinear terms (dashed lines), and 
the trajectories obtained from the improved version of ODEs (continuous lines). In 
particular, the left and the right panels of Fig. 3 show the phases space trajectories for 
the Josephson and self-trapping regimes, respectively, for single component AJJs. For 
two components AJJs, in Fig. 5 we have plotted the phases space trajectories for the 
Josephson regime, and in Fig. 7 we have plotted the phases space trajectories when the 
system is self-trapped. From Figs. 3, 5, 7 we can see that the trajectories predicted 
when the ODEs Eqs. (19) and (29) are solved in the absence of the usually-neglected 
nonlinear terms are sufficiently close to those predicted when these ODEs are solved in 
the presence of the aforementioned terms. Then, the dynamical evolution predicted by 
the standard ODEs reveals to have a good degree of reliability. 

5. Conclusions 

We have analyzed atomic Josephson junctions for a single Bose gas and for binary 
mixtures of bosons in a double-well potential along the axial direction and a strong 
harmonic confinement in the transverse directions. We have shown that for both the 
cases the Hamiltonian belongs to the extended Bose-Hubbard model and besides the 
density-density interaction it contains the pair hopping and collisionally induced hopping 
terms. These terms derive from the overlaps between wave functions localized in different 
potential wells. We started from these Hamiltonian models and established connections 
with spin Hamiltonians. Proceeding from these, we have discussed the possibility to 
discriminate, under certain conditions, different dynamical regimes sustained by the 
bosonic junctions. From the mean field analysis of the equations of motion for the 
single-particle operators involved in the extended Bose-Hubbard Hamiltonians, we have 
obtained the ordinary differential equations that control the macroscopic dynamics 
of the atomic Josephson junctions. Within the analysis of the atomic Josephson 
junctions macroscopic dynamics we have plotted the phase-plane portraits of the 
dynamical variables (fractional imbalance-relative phase) showing that the inclusion of 
the aforementioned collisionally induced hopping and pair hopping terms are crucial to 
get good agreement between the dynamics of the Josephson model described by ordinary 
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differential equations and the one of the time dependent Gross-Pitaevskii equations, 
especially when the atom-atom interaction is strong. 

Finally, it is important to remark that the obtained results are of general validity 
also for more confining (e.g. not saturating to zero at large distances) double-well po- 
tentials. Nevertheless, it is possible to design a model of pair hopping and collisionally 
induced hopping for bosonic atoms that is physically meaningful when optical lattices 
play the role of confining potentials. Physical effects related to pair hopping and colli- 
sionally induced hopping should be observable in generalizations of current experiments 
to detect the superfluid and insulating phases [33]. 

This work has been partially supported by Fondazione CARIPARO through the 
Project 2006: "Guided solitons in matter waves and optical waves with normal and 
anomalous dispersion". G. M. thanks A. B. Kuklov and B. V. Svistunov for useful 
comments. 
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